MODULE menu_contr
  !-----------------------------------------------------------
  ! KRUUNUNHAKA basis set tool kit - contraction optimization menu
  !-----------------------------------------------------------
  ! See CHANGELOG
  !-----------------------------------------------------------
  ! Copyright (C) 2003-2008 Pekka Manninen, 2010 Jussi Lehtola
  !
  ! This program is distributed under the terms of the GNU General
  ! Public License as published by the Free Software Foundation;
  ! either version 2 of the License, or (at your option) any later version.
  !
  ! This program is distributed in the hope that it will be useful,
  ! but WITHOUT ANY WARRANTY; without even the implied warranty of
  ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. 
  !
  ! Please note that this program as well as its 
  ! author must be properly cited whenever the program or some parts 
  ! originated on it are employed. 
  !
  !-----------------------------------------------------------

  USE definitions
  USE textutil
  USE menu

  USE optimize_contractions

  IMPLICIT NONE

CONTAINS

  SUBROUTINE contr_menu(X_l,C_l,nbas)
    REAL(KIND=prec), DIMENSION(:,:), ALLOCATABLE, INTENT(INOUT) :: X_l
    REAL(KIND=prec), DIMENSION(:,:,:), ALLOCATABLE, INTENT(INOUT) :: C_l
    INTEGER, DIMENSION(MAXL,2), INTENT(INOUT) :: nbas

    ! Contraction menu
    CHARACTER(LEN=TERMWIDTH) :: contrmenu_tit
    CHARACTER(LEN=TERMWIDTH), DIMENSION(10) :: contrmenu
    CHARACTER(LEN=10) :: contrsel

    ! The same, when there are no primitives to use
    CHARACTER(LEN=TERMWIDTH), DIMENSION(2) :: contrmenu_noprim
    CHARACTER(LEN=2) :: contrsel_noprim

    ! Selection
    INTEGER :: sel

    ! Current values
    CHARACTER(LEN=TERMWIDTH), DIMENSION(3) :: currentvals
    CHARACTER(LEN=TERMWIDTH), DIMENSION(1) :: currentvals_noprim

    ! Helper array
    CHARACTER(LEN=TERMWIDTH) :: tmpstr

    ! Angular momentum to optimize for
    INTEGER, SAVE :: l=1
    ! Exponent range
    REAL(KIND=prec), SAVE :: low=1e-2_prec
    REAL(KIND=prec), SAVE :: high=1e3_prec
    ! Measure of goodness to use
    INTEGER, SAVE :: measure=2
    ! Tolerance to use in determining convergence
    REAL(KIND=prec), SAVE :: conv=1e-3
    ! Use STO to probe?
    LOGICAL, SAVE :: usesto=.FALSE.

    ! Number of contractions to form
    INTEGER :: ncontr

    ! Loop index
    INTEGER :: j

    LOGICAL :: clearok
    LOGICAL :: finished

    contrmenu_tit="Primitive set optimization menu"

    contrmenu(1)='Set angular momentum type'
    contrmenu(2)='Set range to optimize in'
    contrmenu(3)='Set number of contractions'
    contrmenu(4)='Change measure of goodness of fit'
    contrmenu(5)='Change type of probe'
    contrmenu(6)='Set convergence threshold'
    contrmenu(7)='Show current contraction'
    contrmenu(8)='Remove contractions'
    contrmenu(9)='Create optimized contractions'
    contrmenu(10)='Return to main menu'
    contrsel='123456789r'

    contrmenu_noprim(1)='Set angular momentum type'
    contrmenu_noprim(2)='Return to main menu'
    contrsel_noprim='1r'

    ! Default - don't touch number of contractions
    ncontr=nbas(l,2)    

    finished=.FALSE.
    DO WHILE(.not. finished)
       
       ! If there are primitives on the shell..
       IF(nbas(l,1)>0) THEN
          
          ! Form current values
          WRITE(currentvals(1),'(A,I3,A,A1,A,I3,A)') 'Going to optimize ', ncontr, &
               ' contractions on the ', l_name(l), ' shell which has ', nbas(l,1), ' primitives.'
          
          ! Which measure of goodness is used?
          SELECT CASE(measure)
          CASE(1)
             write(currentvals(2),'(A)') 'Trying to optimize by minimizing maximum deviance'
          CASE(2)
             write(currentvals(2),'(A)') 'Trying to optimize by minimizing rms deviance'
          CASE(3)
             write(currentvals(2),'(A)') 'Trying to optimize by maximizing integral over profile'
          END SELECT

          IF(usesto) THEN
             WRITE(currentvals(3),'(A,E9.2,A,E9.2,A,E9.2,A)') 'Probing with an STO from ',&
                  low, ' to ', high, ' with threshold ',conv,'.'
          ELSE
             WRITE(currentvals(3),'(A,E9.2,A,E9.2,A,E9.2,A)') 'Probing with a GTO from ',&
                  low, ' to ', high, ' with threshold ',conv,'.'
          END IF
          
          
          
          ! Do selection
          sel=make_cvals_selection(contrmenu_tit,currentvals,contrmenu,contrsel,'*')
          
          SELECT CASE(sel)
          CASE(1)
             l=-1
             DO WHILE(l < 0 .OR. l > MAXL)
                l=read_int('Give value of angular momentum (S = 0, P = 1, etc):')
                IF(l < 0 .OR. l > MAXL) THEN
                   WRITE (*,'(A)') 'Invalid value of angular momentum.'
                END IF
             END DO
             ! Convert to Fortran indexing
             l=l+1
             
          CASE(2)
             ! Set limits
             low=1.0
             high=0.0
             DO WHILE(low>high .OR. low<0)
                CALL read_limits('Give limits of scanning range:',low,high)
                IF(low>high .OR. low<0) THEN
                   WRITE (*,'(A)') 'Invalid scanning range.'
                END IF
             END DO
             
          CASE(3)
             ! Set number of contractions
             ncontr=-1
             DO WHILE(ncontr <= 0 .OR. ncontr>nbas(l,1))
                WRITE (tmpstr,'(A,I3,A)') 'Give wanted amount of contractions [1 -', nbas(l,1), ']'
                ncontr=read_int(tmpstr)
                IF(ncontr <= 0 .OR. ncontr>nbas(l,1)) THEN
                   WRITE (*,'(A)') 'Invalid number of contractions.'
                END IF
             END DO
             
          CASE(4)
             ! Change measure of goodness
             measure=0
             DO WHILE(measure <1 .OR. measure>3)
                measure=read_int('1. Minimize deviance, 2. Minimize RMS deviance, 3. Maximize integral?')
                IF(measure <1 .OR. measure>3) THEN
                   WRITE (*,'(A)') 'Invalid measure.'
                END IF
             END DO
             
          CASE(5)
             ! Change probe
             usesto=.NOT. usesto

          CASE(6)
             ! Set convergence threshold
             conv=-1
             DO WHILE(conv<=0)
                conv=read_double('Give value for convergence threshold')
                IF(conv<=0) THEN
                   WRITE (*,'(A)') 'Invalid value for convergence threshold.'
                END IF
             END DO
             
          CASE(7)
             ! Show current primitives
             
             WRITE(*,'(A,A,A)') '$ ', l_name(l), '-TYPE FUNCTIONS'
             WRITE(*,'(2I5)') nbas(l,1), nbas(l,2)
             DO j = 1, nbas(l,1)
                WRITE(*,'(E12.5,A,20F8.4)') X_l(l,j), '  ', C_l(l,j,1:nbas(l,2))
             END DO
             WRITE (*,*) ''
             
          CASE(8)
             ! Remove contractions
             
             clearok=.TRUE.
             
             ! Check if shell already has something
             IF(nbas(l,2) /= nbas(l,1)) THEN
                WRITE (tmpstr,'(A,A,A,I3,A)') 'The ', l_name(l), ' shell has ', nbas(l,2), ' contractions.'
                IF(make_yn(tmpstr)) THEN
                   ! OK to proceed
                   clearok=.TRUE.
                ELSE
                   ! Not ok.
                   clearok=.FALSE.
                END IF
             END IF
             
             IF(clearok) THEN 
                ! Make dummy contraction matrix
                C_l(l,:,:)=0.0
                DO j=1,nbas(l,1)
                   C_l(l,j,j)=1.0
                END DO
                nbas(l,2)=nbas(l,1)
             END IF
             
          CASE(9)
             ! Perform optimization
             
             ! Get affirmation
             WRITE (tmpstr,'(A,A,A,I3,A)') 'The ', l_name(l), ' shell has currently ', nbas(l,2), ' contractions.'
             IF(make_yn(tmpstr)) THEN
                ! OK to proceed
                clearok=.TRUE.
             ELSE
                ! Not ok.
                clearok=.FALSE.
             END IF

             ! The same number of basis functions requested as there are exponents
             IF(ncontr == nbas(l,1)) THEN
                WRITE (tmpstr,'(A,I3,A,I3,A)') 'Really optimize ', ncontr, ' contractions out of ', nbas(l,1), ' exponents?'
                IF(make_yn(tmpstr)) THEN
                   clearok=.TRUE.
                ELSE
                   clearok=.FALSE.
                END IF
             END IF
             
             IF(clearok) THEN 
                ! OK to proceed with optimization

                ! Sanity check wanted number of contractions
                IF(ncontr > nbas(l,1)) THEN
                   ! We actually have checked that it can't be higher than nbas(l,1)
                   ! So we can just make a dummy contraction matrix
                   WRITE (*,'(/A,I3,A,I3,A/)') 'Requested ', nbas(l,2),' contractions, while there are ',&
                        nbas(l,1), ' exponents. Forming dummy matrix.'

                   C_l(l,:,:)=0.0
                   DO j=1,nbas(l,1)
                      C_l(l,j,j)=1.0
                   END DO
                   nbas(l,2)=nbas(l,1)

                ELSE
                   ! Set amount of contractions
                   nbas(l,2)=ncontr
                   
                   !! TODO
                   ! Run optimization
                   CALL opt_contr(l, low, high, X_l(l,1:nbas(l,1)), C_l(l,1:nbas(l,1),1:nbas(l,2)), usesto, measure, conv)
!                  CALL opt_contr_simann(l, low, high, X_l(l,1:nbas(l,1)), C_l(l,1:nbas(l,1),1:nbas(l,2)), usesto,&
!                       measure,5000)
!                   CALL opt_contr_lottery(l, low, high, X_l(l,1:nbas(l,1)), C_l(l,1:nbas(l,1),1:nbas(l,2)), usesto,&
!                        measure, 5000)
                END IF
             END IF
             
          CASE(10)
             ! Return to main menu
             finished=.TRUE.
          END SELECT
       ELSE
          ! No primitives, cannot perform contractions.

          WRITE(currentvals_noprim(1),'(A,A,A)') 'The ', l_name(l), ' shell has no primitives.'

          ! Do selection
          sel=make_cvals_selection(contrmenu_tit,currentvals_noprim,contrmenu_noprim,contrsel_noprim,'*')
          
          SELECT CASE(sel)
          CASE(1)
             l=-1
             DO WHILE(l < 0 .OR. l > MAXL)
                WRITE (*,'(A)',advance='no') 'Give value of angular momentum (S = 0, P = 1, etc): '
                READ (*,*) l
                WRITE (*,*) ''
             END DO
             ! Convert to Fortran indexing
             l=l+1
             
          CASE(2)
             ! Return to main menu
             finished=.TRUE.
          END SELECT
       END IF
    END DO

  END SUBROUTINE contr_menu

END MODULE menu_contr
